Forcing function control of Faraday wave instabilities in viscous shallow fluids 
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We investigate the relationship between the linear surface wave instabilities of a shallow viscous 
fluid layer and the shape of the periodic, parametric-forcing function (describing the vertical ac- 
celeration of the fluid container) that excites them. We find numerically that the envelope of the 
resonance tongues can only develop multiple minima when the forcing function has more than two 
local extrema per cycle. With this insight, we construct a multi-frequency forcing function that 
generates at onset a non-trivial harmonic instability which is distinct from a subharmonic response 
to any of its frequency components. We measure the corresponding surface patterns experimentally 
and verify that small changes in the forcing waveform cause a transition, through a bicritical point, 
from the predicted harmonic short-wavelength pattern to a much larger standard subharmonic pat- 
tern. Using a formulation valid in the lubrication regime (thin viscous fluid layer) and a WKB 
method to find its analytic solutions, we explore the origin of the observed relation between the 
forcing function shape and the resonance tongue structure. In particular, we show that for square 
and triangular forcing functions the envelope of these tongues has only one minimum, as in the 
usual sinusoidal case. 

PACS numbers: 47.35.+i,47.20.-k,47.54.+r 



I. INTRODUCTION 

In the Faraday system, an incompressible fluid is oscil- 
lated vertically in a container with a free upper surface, 
generating standing surface waves which provide an ex- 
cellent system for the study of pattern formation 0, . 
Through an appropriate choice of experimental parame- 
ters, many of the regular patterns that are possible in two 
dimensions, such as stripes, squares and hexagons, can be 
obtained. In addition, targets, spirals, superlattices and 
quasipatterns lacking strict translational periodicity have 
also been observed [3, H IE IS • 

One of the advantages of the Faraday experiment, 
when compared to other pattern-forming systems such 
as convection or chemical reactions, is the great amount 
of control over the energy feeding mechanism that can 
be achieved by changing the periodic vertical accelera- 
tion of the fluid container. Even by forcing the system 
with different combinations of only two frequencies, sev- 
eral distinct patterns can be achieved. Hexagonal and 
rhomboid patterns, together with various quasipatterns 
have been obtained experimentally in [3, HJ, by vary- 
ing the amplitudes and the phase difference between both 
components. Superlattice patterns ITU Il2j . triangular 
patterns and localized structures Il4l have also been 
observed using two- frequency forcings |15|. 

From a theoretical perspective, a combination of tools 
must be used to understand and predict the pattern se- 
lection. While its characteristic wavelength can be ob- 
tained through a linear instability calculation, the two- 
dimensional structure is determined by the nonlinear in- 
teraction between modes [IE E3, EE E3, E3, H3, 123, 
l24j . At the linear level, the simplest cases occur when a 
deep fluid layer of low viscosity is oscillated with a sinu- 
soidal forcing, i.e. proportional to sm(ujt). In these situ- 



ations, the frequency of the main (largest in amplitude) 
component of the resulting surface wave oscillations will 
be cu/2 (referred to hereafter as the first -or fundamental- 
subharmonic response). In other cases, two mechanisms 
for selecting the main frequency responses that are differ- 
ent from the first subharmonic one have been identified. 

The first mechanism occurs when two or more fre- 
quency components are introduced in the forcing. In 
these cases, each component will tend to excite its own 
corresponding first subharmonic mode. Their relative 
amplitudes will determine which of these responses has 
the lowest global forcing strength threshold, thus becom- 
ing the instability that is observed at onset. The second 
mechanism can only arise in the high viscosity regime. If 
the fluid layer is shallow enough, even a single component 
forcing with low enough frequency can excite an instabil- 
ity different from the first subharmonic one. As the vis- 
cous boundary layer reaches the bottom of the fluid con- 
tainer, the threshold of the lowest unstable modes rises, 
allowing others with higher main frequency components 
(and, therefore, shorter surface wavelengths) to become 
unstable at onset p5H2f|. 

In a numerical and experimental study, it was shown in 
[27) that a transition between two patterns with different 
linearly unstable wavelengths can be obtained in various 
fluid regimes by changing the relative amplitudes of a 
two-frequency forcing function. This transition occurs 
through a bicritical point, where both modes are simul- 
taneously neutrally stable. In spite of these results, only 
a limited understanding of the effects of both a multi- 
frequency forcing and a high viscosity regime has been 
achieved. Furthermore, little is known about the pat- 
terns expected for more complicated forcing functions not 
described by a few frequency components. This can be 
attributed to the essentially infinite number of degrees 
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of freedom that are needed to parametrize an arbitrary 
forcing function, which renders a systematic exploration 
of the parameter space impossible. 

In this paper, we consider a different and novel ap- 
proach. Instead of exploring a large parameter space with 
various forcing frequency components, we seek to iden- 
tify which characteristics of the periodic forcing function 
affect the surface patterns and how. By performing a nu- 
merical linear stability calculation in various test systems 
of shallow viscous fluid layers, we will first identify a sim- 
ple qualitative relation between the shape of the forcing 
function and the resonance tongue structure (that de- 
scribes the stability thresholds). Using this relation, we 
will construct a forcing function with a non-trivial critical 
instability at onset, having a main frequency component 
which does not correspond to the fundamental subhar- 
monic (or even the fundamental harmonic) response to 
any of its forcing frequencies. We will then present exper- 
imental results showing the surface pattern generated by 
this instability. Finally, in the lubrication limit of a thin 
viscous fluid layer, we will illustrate analytically the ori- 
gin of the observed relation between the forcing function 
and the stability thresholds. We will follow the method 
introduced by Cerda and Tirapegui [28l I29I ] that derives 
a Mathieu equation to describe this regime and uses a 
WKB approximation [3(3, to solve it for single fre- 
quency forcing. By extending these calculations to arbi- 
trary forcing functions we will develop an intuitive under- 
standing of the relation between the shape of the forcing 
function and the structure of the resonance tongues. In 
particular, we will show that only forcing functions with 
more than two local extrema per cycle are expected to 
allow bicritical points involving non-contiguous tongues. 

The paper is organized as follows. In Section [D] we 
review the standard formulation of the Faraday wave lin- 
ear stability analysis. We introduce in Section ITTT1 a one- 
parameter family of forcing functions to illustrate numer- 
ically the relation between the shape of each member of 
the family and the structure of its corresponding neutral 
stability diagram. Section IIVI presents an experimental 
study that uses these forcing functions, displaying a pre- 
viously unobserved transition between two surface pat- 
terns with very different characteristic wavelengths. In 
Section|V| we show an approximate analytical relation be- 
tween the forcing and the instability response that illumi- 
nates our approach. Finally, Section IVll briefly discusses 
our results and presents our conclusions. 



II. BACKGROUND 

We study the linear stability of the free surface of an 
incompressible Newtonian fluid layer of depth h, density 
p, kinematic viscosity v and surface tension <r. The fluid 
is oscillated vertically with acceleration /(cut), where 10 is 
the fundamental frequency of oscillation and t is the time. 
We will summarize here the derivation of the equations 
describing this system by following the presentation in 



M 

Using the incompressibility condition to eliminate the 
pressure in the linearized Navier-Stokes equation we ob- 
tain 

{d t - vV 2 )V 2 u z = 0, (1) 

where u z (x, y, z, t) is the vertical component of the fluid 
velocity. In an idealized laterally infinite container, the 
horizontal eigenfunctions are given by e ±lk ' r , with r = 
(x,y) and k — (k x ,k y ). For each surface wavenumbcr 
k = \k\, equation thus becomes 

[dt ~ v{d„ - fc 2 )] (d zz - k 2 )v k = 0, (2) 

where Vk(z,t) describes the z-dependence of u z associ- 
ated with the mode k. In the oscillating reference frame 
with z = at the flat fluid surface, the boundary condi- 
tions on the bottom of the container are given by 

Vk = and d z v k = 0, at z = —h. (3) 

At the fluid surface, the vertical position of the free 
boundary z — £fc(i)e lfe ' r associated to every mode k is 
advected by the fluid motion. To linear order in the sur- 
face deformation this kinematic boundary condition is 

^ = v k at z = 0. (4) 

Additional boundary conditions are imposed at the sur- 
face by finding the total balance of forces tangential and 
normal to the interface. From this we obtain 

(d zz + k 2 ) v k = 0, (5) 
[dt - v[d zz - k 2 ) + 2vk 2 } 8 z v k = 

[ fl (i + r/(wt)) + p -k 2 \k 2 i k , (6) 

at z = 0, 

where g is the gravitational acceleration and / is a non- 
dimensional function defined to have max(|/(u;t)|) = 1. 
Therefore, T corresponds to the maximum acceleration 
of the forcing function, expressed in units of g. 

Equation © and boundary conditions © through © 
fully describe the dynamics of the system. Instead of 
integrating them directly, our numerical analysis will fo- 
cus on finding the stability threshold T c (k) given by the 
critical value of T at which the wavenumber k becomes 
unstable. 



III. NUMERICAL STUDY 
A. Method 

We are interested in finding numerically the neutral 
stability curve T c (fc) for various forcing functions. With 
this objective, we have extended the stability analysis 
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method of Kumar and Tuckerman [27II32] ] to forcing func- 
tions with an arbitrary number of frequency components. 
In broad terms, this method consists first in expanding 
Vk and in a Floquet form 

Vk = e (M+#)< J2 Wj {z)e ijut + c.c. (7) 

3 

a- = e^+^t Qe l]UJt + c.c. (8) 

3 

Here, p, + is the Floquet exponent, where we can set 
the growth rate p, to to obtain marginal stability curves 
with harmonic (<f> — 0) and subharmonic (</> = lu/2) 
temporal responses. Equations ©-© are then used to 
rewrite © in the form 

AnCn = T [/C]„ , (9) 

where A n is an algebraic function of the system param- 
eters, which does not depend on f(u>t), and [/C]„ is the 
n-th Fourier component of 

f(u;t)J2Cje ijut . (10) 

3 

By introducing the explicit form of /(cot), equation @ 
can be expressed as an eigenvalue problem for the forcing 
amplitude T which can then be solved through standard 
numerical techniques. In order to extend the method 
to cases beyond the two-frequency forcing computed in 
|27l l32| , we implemented this algorithm in Mathematica 
[33| and used the program's symbolic algebra capabili- 
ties to automatically compute [/C]„ for any given f(ujt). 
With this implementation, which is analogous to that 
presented in [34(, we are able to obtain efficiently the 
neutral stability curves for any desired forcing function, 
regardless of its frequency content. 

B. Results 

We restrict our study to shallow viscous fluid layers. 
Since the specific value of the fluid constants within this 
regime does not change our qualitative results or analysis, 
we will further reduce the size of the parameter space by 
considering throughout the paper only one set of fluid 
constants. These are given by a density p — 0.95 g/cm 3 , a 
surface tension a = 20dyn/cm and a viscosity v — 46 cS. 
Additionally, we will use in this section and in Section 
HVI a fluid depth h — 0.3 cm and an oscillation frequency 
w = 2tt(10Hz). 

By using the numerical techniques described above, we 
explored the structure of the marginal stability curves 
T c (k) for many different f(ujt) including various piece- 
wise constant, piecewise linear, delta-like and multi- 
frequency functions. While a precise characterization of 
how the features of f{tot) correlate to those of T c (k) re- 
mains to be achieved, one of the salient qualitative rela- 
tions that we observed for all tested functions is a con- 
nection between the extrema of f{ut) and the envelope 




FIG. 1: (Color online) Shape of the forcing functions (left) de- 
fined in pip and their corresponding neutral stability curves 
(right) for (a) p = -2, (b) p = -0.3, (c) p = 0.5, (d) p = 1, 
and parameters p — 0.95 g/cm 3 , a = 20dyn/cm, v — 46 cS, 
uj — 27r(10Hz) and h = 0.3 cm. F is in units of g and k 
in cm' 1 . The resonance tongues labeled H and SH show 
regions with harmonic or subharmonic linear instabilities, re- 
spectively. Note how their envelopes (dashed lines) change 
with p. 

of r c (fc) that will be described below. We will illustrate 
it here for a specific family of forcing functions, which is 
the same as that in the experiments of Section IIVI 

Consider the following set of forcing functions 
parametrized by p 

f p (ujt) = Af [2.5 cos(wi) + 3 P cos(3wi) - 5 P cos(5w£)] , 

(11) 

where to is the fundamental frequency of oscillation and 
j\f is a normalization constant which is defined so that 
max(|/ p (w£)|) = 1. The specific form of pip is an ar- 
bitrary choice which is not important for the qualitative 
behavior that we will focus on here. It was obtained 
by searching for a one-parameter family of forcing func- 
tions that simultaneously includes members with a sim- 
ple triangular-like form (p w —2) and others that can 
produce non-trivial surface-wave instabilities in an ex- 
perimentally accessible regime (p ~ 1). 

Figure n displays in the left column f p (tot) for p = 
—2, p = —0.3, p — 0.5 and p = 1. The right col- 
umn shows the corresponding neutral stability curves 
F c (/c) which present the usual resonance tongue struc- 
ture. The harmonic and subharmonic tongues indicate 
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regions where surface waves become unstable, oscillat- 
ing with a main frequency component that is an in- 
tegral multiple (u),2u),3u>, . . .) or an odd half- multiple 
(w/2, 3w/2, 5w/2, . . .) of the fundamental forcing fre- 
quency, respectively. The tongues at higher fc-values cor- 
respond to instabilities with shorter surface wavelengths 
and higher oscillation frequencies. As p is increased, the 
forcing function changes from a simple rounded triangu- 
lar shape with only two extrema per cycle to shapes with 
richer structure. Correspondingly, the envelope defined 
by the tongue minima (sketched as a dashed line on the 
figure) changes from a simple convex function with a sin- 
gle minimum to a set of convex segments, each with its 
own minimum. 

We have observed a similar relation between the struc- 
ture of the extrema of f(tot) and the concavity of the res- 
onance tongue's envelope for all forcing functions tested 
(triangular, square, multi-frequency, etc.) In particular, 
every f{uit) with only two extrema per cycle resulted in 
an envelope with positive concavity for all k. This rela- 
tion will be one of our main focuses in the remainder of 
this paper. 

It is important to point out that the changes in the 
critical instabilities illustrated in Fig. ^ cannot be ex- 
plained by a simple switch to a different dominant forcing 
frequency in f p (u)t) combined with the first mechanism 
described in the Introduction. Indeed, as p is increased 
to 1 the lowest unstable region becomes the second har- 
monic tongue (with main frequency component equal to 
2lu) which does not correspond to the fundamental har- 
monic or subharmonic responses (with equal or half the 
frequency, respectively) to any of the three frequency 
components of f p {ut): u>, 3lu and 5lo. Furthermore, it 
is apparent that the change in p cannot be characterized 
as mainly reducing the stability threshold of a specific 
tongue, but that it rather affects the aforementioned en- 
velope over the entire range of k studied. 



IV. EXPERIMENTAL RESULTS 

In this section, we present experimental results show- 
ing that the appearance of multiple minima in the en- 
velope of the resonance tongues can generate interesting 
measurable effects. By carefully choosing the form of the 
forcing function, we find a previously unobserved bicriti- 
cal point between two surface patterns with very different 
characteristic wavelengths. 

In our experiments, we use silicone oil with p = 
0.95 g/cm 3 , a = 20dyn/cm and v = 46 cS (Fluka Sili- 
cone Oil AR 20), which are the same fluid parameters 
as in Section ITTll A 0.3 cm deep layer of this silicone oil 
is contained in a cylindrical cell with a radius of 7.0 cm 
and height of 4.0 cm. The cell has a PVC sidewall, a 
0.8 cm thick glass bottom, and a 0.8 cm thick plexiglass 
top covered with a light diffuser. It is mounted on the 
ram of a 10 x 10 cm 2 linear air bearing, which is attached 
to a 180 kg triangular granite slab that floats on an air 



table to minimize horizontal oscillations. A shaker (VTS 
VG100) is suspended by springs from the air table sup- 
ports. Two 50 cm long cylindrical aluminum tubes, each 
with an inner and outer diameter of 0.48 cm and 0.95 cm, 
respectively, connect the shaker to the ram. An ampli- 
fier (Crown CE2000) drives the shaker with a computer- 
generated forcing function. The amplitudes and phases of 
the desired Fourier components of the acceleration signal 
are measured by an accelerometer (PCB Model 353B68) 
and used as feedback to control the driving. The root- 
mean-square difference between the measured and target 
forcing functions is less than 1% while the variation in the 
amplitudes of the driven components is less than 0.01%. 
Since viscosity and surface tension are both sensitive to 
temperature changes, the experiments are conducted in a 
closed transparent box maintained at a constant temper- 
ature (±0.005°C). To visualize the waves, parallel light 
is projected through the cell bottom. The curved fluid 
surface refracts the light, which then falls on the diffuser 
producing a representation of the pattern. A CCD cam- 
era synchronized with the forcing function acquires the 
images. 

Our specific choice of forcing function was determined 
by searching for an experimentally achievable set of pa- 
rameters having a linear instability at onset with a re- 
sponse far from the usual subharmonic one. This objec- 
tive is not easily achieved despite the fact that our nu- 
merical exploration established that many forcing func- 
tions generate resonance tongues with a multiple minima 
envelope. Indeed, for the fluid parameters used in our ex- 
periments, we found numerically that a tongue belonging 
to the second or higher (in order of increasing k) enve- 
lope minimum can be excited at onset only for very low 
values of h or uj. However, the range of these two quan- 
tities is limited by our experimental apparatus. For very 
shallow fluid layers (h < 0.1 cm), spurious effects can af- 
fect the patterns: surface waves may contact the bottom 
of the container and a small tilt, variation in the bottom 
profile, or wetting at the wall can lead to large changes 
in the relative fluid depth Ah/h. Additionally, as h and 
u> are reduced, the critical acceleration r c increases. Be- 
cause the maximum acceleration and amplitude (oc uj~ 2 ) 
of the apparatus are limited, much of this low w/large 
r regime is inaccessible. By testing numerically various 
forcing functions, we were able to construct f p (u>t) with 
p w 1 , as defined in (| 1 1 1> , which has a global minimum in 
the part of the envelope that does not contain the first 
subharmonic tongue (see Fig. [fll), and which is experi- 
mentally accessible. 

Figure|5]displays the neutral stability curves computed 
numerically for the experimental parameters specified 
above, using to = 27r(10Hz) and a forcing f p (u>t) with 
p = 0.9, p = 1.0 and p = 1.1. The figure shows a 
very small change in the forcing function (see left panels) 
producing a large jump in the critical wavenumber. For 
p = 0.9 (top), the first subharmonic tongue (with main 
frequency component at ui/2) will be excited at onset. 
Numerically, we compute a critical forcing = 3.35 
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FIG. 3: Experimental pictures (negative-images) of the sur- 
face patterns appearing at onset for p = 0.9 (left) and p = 1.1 
(right), corresponding to the top and bottom forcing functions 
in Fig. |21 Note that, for this small variation in the forcing, a 
dramatic change in the pattern is observed. The size of each 
image is 8.22 cm x 8.22 cm, which captures the central region 
of the container. 
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FIG. 2: (Color online) Shape of the forcing functions (left) 
and their corresponding neutral stability curves (right) for 
p = 0.9 (top), p — 1.0 (center), p — 1.1 (bottom), and the 
same parameters and units as in Fig. Q For small changes 
in f(uit), the instability at onset (occurring at critical forc- 
ings indicated by the dashed lines) switches from the first 
subharmonic (SH) resonance tongue to the second harmonic 
(H) one. Each asterisk indicates the critical forcing F c and 
wavenumber k c measured experimentally. 



and a critical wavenumber fcf H = 2.07. At p = 1.0 (cen- 
ter), the system is close to a bicritical point, where the 
first subharmonic and second harmonic tongues become 
simultaneously unstable at onset. The corresponding 
critical values are Tf 1 = 3.86, fcf H = 2.05 and Tf = 3.87, 
= 7.64, respectively. Finally, for p = 1.1 (bottom) 
the second harmonic tongue becomes the instability at 
onset, with = 4.10 and kf = 7.59. We refer to it as 
the second harmonic one since it oscillates with a main 
frequency component at 2tu, and is therefore the second 
harmonic tongue in order of growing k (the first being 
above the plotted T-range, between the two subharmonic 
tongues displayed). It is a non-trivial critical instabil- 
ity, which cannot be easily explained by the mechanisms 
described in the Introduction, because it does not corre- 
spond to the first harmonic or subharmonic responses to 
any of the forcing frequency components (cj, 3lu and 5w). 
Instead, it is related to the second local minimum of the 
envelope of the resonance tongues (see Fig. ^1) . 

Experimentally, we observe the transition between 
these two linearly unstable regimes by using the same 
f p (uit) forcing function. Figure [21 shows images of the 
surface patterns for p — 0.9 (left) and p = 1.1 (right). 
As predicted by our numerical calculations, their char- 
acteristic length scale changes dramatically, in spite of 
the small variation in f p {ojt). For p — 0.9, we obtain a 
pattern of large hexagons at a critical forcing T c — 3.72g, 
with a characteristic size of 3.5 cm which corresponds to 



the critical wavenumber k c — 1.77. When compared to 
the numerical results, T c is within 10% and k c within 
15% of the predicted values. Given the pattern defor- 
mation that is observed towards the image borders due 
to the small aspect ratio (the size of the container is 
only about twice the surface wavelength), these discrep- 
ancies are not significant. For p = 1.1, a pattern of small 
hexagons appears at T c — 4.27, with a characteristic size 
of 0.9 cm, which implies k c = 6.96. These measurements 
are within 4% (for T c ) and 9% (for k c ) of the numeri- 
cal predictions. We have also verified in our experiments 
that, with respect to the fundamental forcing frequency, 
the oscillations of the large pattern are subharmonic and 
those of the small one are harmonic. Finally, at p = 1.0 
(image not shown), we observe that the system gener- 
ates small hexagons which are practically indistinguish- 
able from those at p — 1.1, with r c = 4.0 and k c = 6.98. 



For 0.92 < p < 0.95, we find in our experiments a bi- 
critical region where a complicated mixed mode surface 
pattern appears. These kind of patterns can arise from 
the nonlinear interactions of two or more linear insta- 
bilities 

[H SJ El S H|. They are often obtained by 
introducing frequency components in the forcing function 
with simple linear responses that interact in the horizon- 
tal plane to produce new structures. In contrast, in the 
current situation the changes in the tongue envelope se- 
lects linear instabilities that are not directly connected 
to the forcing components, and therefore the patterns 
generated through this mechanism could potentially be 
different. Unfortunately, in our current experiment the 
mixed surface patterns include complicated interactions 
with the side walls due to the small size of the container. 
Their proper analysis will therefore require a much larger 
aspect ratio and is left for future work. 
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FIG. 4: (Color online) Plot of functions Bi (solid line) 
and B2 (dashed line), as denned by equations lib! and 1171 , 
Both functions diverge for y — > with Bi(y) <x y~ 3 and 
B2(y) tx For y — > 00 they approach their corresponding 
asymptotic limits Bi(y) — > 1 and B2(y) — ► 3/2. 



V. ANALYTICAL CALCULATIONS 
A. The lubrication approximation 

We are interested in exploring analytically the origin 
of the relation observed in Section ITTT1 between fiuit) and 
the envelope of r c (fe). To proceed, we will focus on sys- 
tems in the lubrication regime, where the ratio between 
the dtv and the V 2 ?7 term of the Navier-Stokes equation 
is small. This ratio is of order (l/S) 2 , where I is the dis- 
tance that the fluid motion penetrates the surface and 8 is 
the characteristic size of the boundary layer [2^, I29I [35J . 
Since / can be estimated by either 1/k (if kh ^> 1) or 
h (if kh < 1), and S is proportional to yV/u;, it follows 
that a system is in the lubrication regime if it consists of 
a shallow enough fluid layer with high enough viscosity 
and a low enough oscillation frequency. 

We use a simplified anal ytic description, introduced 
by Cerda and Tirapegui in 28, 29] for fluids under the 
lubrication approximation, in which a damped Mathieu 
equation involving only the motion of the free fluid sur- 
face is obtained. This equation is found by first deriving 
an exact non-local (in time) relation for the linear evo- 
lution of the surface, which is a formulation analogous 
to that developed in [3(|. By imposing a short-memory 
to the system due to its fast dissipation rate, the non- 
local dependence is then neglected. The resulting Math- 
ieu equation reads 



& + 27fcC + Wfe [1 + r fc /(wt)] & - 0, 



(12) 



where the dots represent derivatives with respect to time 
and 



7 fe = vk 2 B 1 (kh)B 2 {kh) 
Q 2 = k [g + ak 2 /p] B 2 (kh) 

r 5 



r fc = 



g + ok 2 1 p 



(13) 
(14) 

(15) 



Here, B\{kh) and B 2 (kh) are explicit non-dimensional 
functions given by 



Bi (y) 



cosh(2t>)+2t( +1 
sinh(2iy)-2y 



, . _ 3cosh 2 (a)[ S inh(2y)-2y-4 a 3 /3]+y 2 [sinh(2y)-2 a ] 
a2 \y) ~ [sinhr2i,1-2j/] 2 



(16) 
(17) 



[sinh(2y)- 

Figure0|shows that Bi(y) and B 2 (y) have a simple struc- 
ture despite their complicated algebraic expressions. As 
y approaches 0, both functions diverge with B\{y) oc y" 3 
and B2{y) oc y~ l . For large values of y, B 1 (y) and B 2 {y) 
quickly converge to their asymptotic limits of 1 and 3/2, 
respectively. 

The critical forcing strength r c can be found for every 
k by considering solutions of (|12|l that follow the Floquet 
form 



&(t + 27r/w) 



= e^~ 



■%<!>) 



(18) 



and demanding that the growth rate after every period 
satisfies p = 0. 



B. The WKB approximation 

We will follow here the approach in psl . |29| . which 
uses the well known (in the context of quantum me- 
chanics) Wentzel-Kramer-Brillouin (WKB) approxima- 
tion |30L l3ll | to solve the Mathieu equation. We first cast 
(|12fl into the form of a Schrodinger equation by defining 



x = tut 
9{x) = t k {x/u)e^ x '« 



and 



to obtain 



E 
V(x) 

1 



-r k uj 2 k f(x), 



*"(x) + -=[E- V(x)} = 0, 



(19) 
(20) 

(21) 
(22) 

(23) 



where the double prime represents the second derivative 
with respect to x. The problem of finding the solutions 
of i|12|) that follow the Floquet form l |18|l then becomes 
equivalent to finding the eigenfunctions of 1)23(1 that sat- 
isfy 



tf (a; + 2tt) = e^ + ^+^<J<(x), 



(24) 



where the neutral stability curves are obtained for p = 0. 

In regions where cu 2 /\E~ V(x)\ <§; 1, the WKB approx- 
imation provides explicit solutions for (|23|l which are di- 
vided into two different families. For E < V(x) (as in the 
(cij, &j)-intervals of Fig. |SJ) they are given in their most 
general form by 



*(x) = 



^PJx) 



Aexp 



P(x)dx 



B exp 



P(x)dx 



(25) 
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u 






with 



% ft* 



FIG. 5: Illustration of the intervals of V(x) in which the 
integrals a-, and /3j are computed using expressions 1311 and 
1321 . In the WKB approach described (see text), a matrix Mj 
is defined through expression I3U1 for each (ctj, Qj + i)-interval. 



and for E > V(x) (intervals (bj,aj+\) in Fig. EJ), by 
1 



C sin 



Z? cos 



P(x)dx 



(26) 



Here, P(;e) = y/\E — and the complex constants 

A, -B, C and D are obtained by imposing the boundary 
conditions in each segment. The solution for a given V(x) 
over the full x domain is found by matching adjacent 
segments of ^(x) at the points xq where V(xq) = E. 
For x w xq, however, expressions <|25[1 and i|26[l are not 
valid and, following the WKB method, one must perform 
a matched asymptotic expansion around xq to find the 
correct matching formulas j^H. At the points {&j}j=i..jv 
shown on Fig. [SJ these are given by 



C = 



2A-B 

~7T~ 



and D = 



2A + B 



and at the points {cij}j 



.N 



by 



A C + D 1 o 

A = =— and B = 



V2 



2V2 



(27) 



(28) 



We will now extend the neutral stability calculations 
carried out in [2^, [2j| for V(x) cx cos(x) to arbitrary 
forcing functions. Imagine a periodic function V(x) with 
2N matching points per period as in Fig. [5] Using 127|) 
and (|28|l we can relate the coefficients A,+i and -Bj+i of 
solution l|25|) in an interval (<ij+i, to the coefficients 
Aj and Bj in the previous interval (a,j,bj) (see Fig. |5J|. 
We find 



B H 



where the matrix Mj is defined by 



Mj = 



e aj sin(/3j-) 



-e aj sin(/3j) 
■e~ Q ^ cos(/3j) 



(29) 



(30) 



n, = / P{x)dx 



ft 



P{x)dx. 



(31) 
(32) 



The change in the amplitude of the wave function ^(a:) 
after a full period is therefore given by the product M = 
M^Mjv_j . . . Mi. Hence, for solutions with the Floquet 
form, equation (|24|l implies the neutral stability condition 



max(|A+|, |A_|) = e" 



(33) 



where a+ 



A + and A_ are the two eigenvalues of M. An 
equivalent condition can be found by using the fact 
that the trace Tr(M) is real and that the determinant 
Det(M) is equal to 1, together with the standard rela- 
tions Tr(M) = A+ + A_ and Det(M) = A+A_. The 
resulting expression is 



2tt 

Tr(M) = ±2 cosh ( — % 

1 UJ 



(34) 



where the plus or minus signs provide the neutral stabil- 
ity boundaries for harmonic or subharmonic resonances, 
respectively. 

Note that for some values of k and T it is also possible 
to have E > V(x) or E < V(x) for all x, and therefore 
no intersections between V(x) and E. In these situations 
the matrix M cannot be computed and our current im- 
plementation breaks down. However, the WKB method 
is still valid and it has been shown in |28l |29| that these 
cases never lead to instabilities. In our computation of 
the neutral stability curves we can therefore assume that 
there is at least one a and one (3 region per cycle. 



C. Validity of the approximation. 

We will investigate here the validity conditions for the 
approximation described above. The WKB method is 
based on an expansion in the small quantity oj 2 /\E — 
V(x)\ which can be estimated by |28ll29| 



\E-V(x)\ 



(35) 



This criterion implies that the approximation should be 
valid for systems with (l/S) 4 -C 1, which is a condition 
that must be satisfied in the lubrication regime in which 
we are focusing. Indeed, the lubrication regime requires 
(l/S) 2 <C 1 and therefore, given that (l/S) 4 will be even 
smaller, the WKB approximation must also be valid in 
this regime. Let us estimate S and I for the fluid parame- 
ters used in Sections ITTT1 and llVl For surface waves oscil- 
lating at a frequency , the characteristic size S of the 
viscous boundary layer is of order ^v/Q.^ psl I2H l35| . 



8 






2 4 6 S 10 2 4 6 S 




FIG. 6: (Color online) Neutral stability curves for a forcing 
function f p (u}i) with (a) p — —2, (d) p = 1 (labeled as in 
Fig. and parameters p — 0.95 g/cm 3 , a = 20dyn/cm, v = 
46 cS, ui — 27r(3.5Hz) and h = 0.1cm. T is in units of g 
and k in cm - . The exact numerical computations (top) are 
compared to the WKB approximation (bottom) . The shape of 
the harmonic (H) and subharmonic (SH) resonance tongues is 
essentially identical for p — — 2 and has similar characteristics 
for p = 1. In both cases, the tongues that would become 
unstable under a forcing of F* = 18 (dashed line) coincide. 



-1! 




FIG. 7: (Color online) Integration regions for the WKB cal- 
culations pertinent to Figs. |5] and |5] The rescaled V(x) = 
V(x)/(TkU>k) curves correspond to a forcing f p (x) with (a) 
p = —2 and (d) p = 1 (labeled as in Fig. 0. The values 
of Ek = E I (r k&k) are displayed for a forcing strength of 
T* = 18 (see Fig. |HJ at k = 4, k = 6 and k = 8. Note that 
the integration zone af is not present for Ee and Eg in (d). 



Since the response frequency of the dominant surface 
waves is typically of the same order as the forcing fre- 
quency, we have that 8 ~ ^/0. 46/10 w 0.2 cm. On the 
other hand, the distance I that the motion of the sur- 
face penetrates the fluid can be estimated by the smallest 
value between h = 0.3 cm and l/k. In the region of k con- 
sidered (see Fig. I is therefore larger than ~ 0.1cm. 
Hence, for these parameters we have that 1/8 is of order 
1, which implies that the WKB method does not provide 
a good approximation. 

In order to be able to use a WKB analysis in our study, 
we will consider in this section a shallower fluid layer with 
h = 0.1cm and a lower oscillation frequency of 3.5 Hz, 
while keeping all other parameters unchanged. For this 
case, we have 8 ~ ^0.46/3.5 w 0.4, and I ~ h = 0.1. We 
thus obtain (1/8) 4 < 10~ 2 , which should imply a good 
WKB approximation. However, this criterion alone does 
not guarantee the accuracy of the resulting neutral sta- 
bility curves. Indeed, for any forcing function there will 
be regions of x where ui 2 /\E — V(x)\ ^S> 1, in which 125|) 
and {25 

are not good approximations. Unfortunately, 
the effect of these regions over the full periodic ^f(x) so- 
lution cannot be easily estimated. This problem becomes 
even harder if V(x) has a complicated shape because in 
such cases no simple approximation can even provide the 
number or size of these regions, which depend on k and 
r. We will therefore validate our analysis by directly 
comparing the WKB results to the numerical solutions 
of the full Navier-Stokes linear stability problem. 

Figure shows the neutral stability curves obtained 
using ui = 2tt (3.5 Hz), h — 0.1cm and the forcing func- 
tion fp(u>t) defined in expression i|ll|) with p = —2 and 
p = 1 (labeled here 'a' and 'd', as in Fig. QJ. The top 



panels show the exact numerical results computed us- 
ing the method described in Section UTTI while the bot- 
tom ones present the approximate WKB solutions. The 
implementation of the WKB algorithm consists in find- 
ing the values r c (fc) for which the trace of M satisfies 
(|34|l . where M is obtained by multiplying the explicit ex- 
pressions for Mj given in l|30l) . By comparing the top 
and bottom panels, it is apparent that the WKB curves 
are almost indistinguishable from the exact results in the 
p = — 2 case. For p = 1, the WKB approximation and 
the exact solution present a similar tongue structure but 
they do not coincide in the exact predicted values for 
the critical stability threshold of each tongue. However, 
the characteristics of their resonance tongue envelopes 
are the same. This is the relevant feature here since it 
is this envelope structure that we will study below using 
the WKB method. 



D. Analysis of the envelopes 

Using the WKB approximation, we are now in a po- 
sition to relate the shape of the forcing function to the 
resonance tongue envelope. For any k and T, the stability 
criterion (|34|) can be computed in terms of 



Q(k,T) = ± 



Tr(M) 



2 cosh (2njk/to) 



(36) 



where Q(k, T) > 1 indicates an instability. If the forc- 
ing function has only two extrema per cycle, there will 
always be at most one a and one (3 integration region, 
as illustrated on Fig. (top) for f p (ojt) with p — — 2 (la- 
beled by an 'a', as in Figs. nand|SJ|. In these cases we 



FIG. 8: (Color online) Plot of Q a (k, T*) for a forcing function 
f P {iot) with p = —2 and a forcing strength F* (see Fig. EJ- The 
regions with Q a > 1 are unstable with harmonic (solid curve) 
or subharmonic (dashed) responses. The dotted envelope is 
computed by discarding the cos(/3 a ) factor in equation 1371 1. 



FIG. 9: (Color online) Plot of Qd(k, T*) for a forcing function 
f p (ujt) with p — 1 and a forcing strength T* (see Fig.EJ- The 
regions with Qd > 1 are unstable with harmonic (solid curve) 
or subharmonic (dashed) responses. At k the definition of 
Qd switches from to , given by Eqs. 138H and 143L 
respectively, since the integration region af is not present for 
k > k (see Fig. 01). 



have M — M\, and (|3"fjj) becomes 



Q a (k,T) = ± 



cosh(cv Q +log2)cos(/3 a ) 
cosh (2tt%/oj) 



(37) 



If we consider the function Q a {k) at constant T, the 
cos(/3 a ) factor will be responsible for oscillations that gen- 
erate an unstable tongue at every excursion that reaches 
Q a > 1. Figure plots Q a at a fixed forcing strength 
r* = 18 g, indicated by the dashed horizontal line on 
Fig. El The dotted lines trace the envelope of Q a , which 
is readily obtained by discarding the cos(/3 a ) factor from 
(|37|l . It exhibits a single maximum on the figure and for 
all other values of T tested, implying that the envelope 
of the resonance tongues must have a single minimum. 

In contrast, forcing functions with multiple extrema 
produce more complicated envelope structures. Figure 
El shows a plot of Q(k, T*) for f p (tot) with p = 1 (la- 
beled here Qd since it corresponds to case 'd' in Figs. 
El and 0) . The oscillation amplitude presents two dis- 
tinct zones of local maxima at k ~ 4 and k ~ 7, which 
are responsible for the two minima that the envelope of 
the resonance tongues displays in Fig. El In general, it 
is easy to see that any resonance tongue envelope with 
multiple minima must be associated with Q(k) functions 
(at fixed r values) which have amplitude envelopes with 
multiple maxima. We will now study how these compli- 
cated amplitude envelopes arise by examining in detail 
the analytical form of Qd- 

The bottom panel of Fig. shows the integration re- 
gions for the p = 1 case. Here, M is given by the prod- 
uct of either three or four matrices, depending on the 
fc-interval considered, since the region is present for 
k < k w 5.7, but not for k > k. In the k < k case it is 
straightforward to compute that 



with 



H<(k) = 
H<(k) = 



cosh(c^ + 2cv^ + af + log 16) 

cosh(27T7fc/o;) 
cosh(af + - ol\ + log 4) 

cosh(27T7fe/ti;) 



and 



C<(k) 
S<{k) 



cos 2 (#) 
-sin 2 (pf) S<(k) 



,2 



(39) 
(40) 

(41) 
.(42) 



In H38|) . we have neglected several additional terms of a 
similar form, but where the argument of the hyperbolic 
cosine contained — af or —a^ contributions. These terms 
turn out to be negligible when compared to H^{k) and 
Hg(k) since a\ and «2 are of the same order, and are 
much larger than a 3 (see Fig. 01) . 

For k > k, M is composed of the product of only three 
matrices and the expressions become simpler. Using an 
equivalent approximation we obtain 



with 



H>(k) 
C>(k) 



cosh(cv 



til 



2a d 2 



log I 



cosh(27T7fc/o;) 
cos 2 (Pi) C>(k) = cos (Pi) 



(43) 

(44) 
(45) 



He ^2 + Hf ( 38 ) 



Figure ITU1 plots the H, C and S functions given above. 
After close examination, one finds that the structure of 
the envelope of Qd(k) is more complicated than that of 
Q a {k) mainly because of the interplay between the oscil- 
lating C and S terms. Indeed, the hyperbolic H terms 
behave similarly to the Q a (k) case, presenting only one 
local maximum, and are therefore not directly related to 
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with 
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FIG. 10: (Color online) Main sinusoidal (top) and hyperbolic 
(bottom) components of the Qd(k,T*) functions displayed in 
Fig. [5] Their combination through equations I38H and I43H 
(for k < k and k > k, respectively) determines the amplitude 
envelope structure observed in Fig. [5] 



the appearance of multiple extrema in the envelope. For 
example, at k w 5 both H^(k) and Hg{k) grow with 
k but the envelope of Qa{k) decreases, mainly because 
of the oscillations of the C^C 2 product. Note that the 
change in the number of integration regions at k is not es- 
sential either for obtaining multiple extrema: the combi- 
nation of the oscillations of the C and S functions are able 
to produce additional extrema even beyond their corre- 
sponding domains. Furthermore, in various tested cases 
with different fluid parameters and forcing functions we 
have found no clear correlation between the changes in 
the number of integration regions and the shape of the 
neutral stability curves. 

We now find analytic expressions that describe the 
envelope of the resonance tongues for any forcing func- 
tion with only two extrema per cycle. The neutral sta- 
bility criterion in these cases is equivalent to setting 
Q a (k, r) = 1 in expression (|37[) . By dropping the oscilla- 
tory factor cos(/3 a ) in (|37[) and using the high dissipation 
of the lubrication regime to neglect the log 2 term (when 
compared to a a which, for the parameters used in this 
section, is evaluated as a a « 27T7fc/u; w 300), we find 
that 



27T7A: 



(46) 



at the envelope. Using the definitions of a a and jk, this 
condition can be rewritten as 



V|l-X + T fc (x)f(x)\dx = 2irjx, (47) 

V(x)>E 

where the integration is carried out over the V(x) > E 
region and the algebraic function x(fc/i) is given by 



x(y) 



^- 2 Bl{y)B 2 {y), 



1 + n 2 y 2 



(48) 



gh 3 



and K2 



gph 2 



(49) 



Equation (|47|l provides an implicit expression for Ffc(x) 
at the envelope. Using this result and the definition in 
l|15p. we find that the shape of the envelope of the res- 
onance tongues under the current approximations is de- 
scribed by the function 



T e {kh)= [1 + K2 k 2 h 2 ] T k ( x (kh)) 



(50) 



Unfortunately, there appears to be no simple way to ex- 
tract the properties of T e (kh) without further specifying 
«i, «2 and f[x). However, we have observed for all tested 
cases that if f(x) has only two extrema per cycle, T e (kh) 
has only one minimum. While the validity of this state- 
ment for all cases is a conjecture that would require a 
proof which is beyond the scope of this paper, we con- 
sider below two simple examples where analytic progress 
can be made. 

For square forcing (where f(uxt) = 1 during half of 
the period and f(u>t) — — 1 during the other half), the 
conjecture can be proved as follows. First, we find the 
solution of (|T7I) 



r fc (x) = 3 X + i. 



(51) 



Then, we substitute this result into equation l|50fl to ob- 
tain an explicit expression for the envelope of the reso- 
nance tongues 



T s e q {kh) = 3 Kl k 6 h 6 B(B 2 + K 2 k 2 h z + 1 



(52) 



While the specific form of T| q depends on the parame- 
ters Ki and K2, its extrema can be readily computed by 
using dkT s e q — 0. We find that they are located at the 
intersection of the functions r(y) — —3 d y [y 3 B 2 (y)B 2 {y)] 
and s(y) — 2k2?//ki. Given that r(y) does not depend on 
any parameters, it can be evaluated numerically without 
loss of generality. We find that it decreases monotoni- 
cally, intersecting the r = axis at y* w 1.479. Using 
this result and the fact that s(y) is a linearly increasing 
function, it is easy to see that T s e q (kh) can have only one 
minimum (which must be located at k < y* /h). 

For triangular forcing, (where f(cot) is a linear function 
that increases during half of the period and a decreases 
during the other half), the analytical calculation becomes 
much harder. The solution for is given by the real root 
of the cubic equation 



(r fc + x - 1) = 9r 2 feX . 



(53) 



It has a more complicated structure than (|51|l . which ren- 
ders the use of the techniques developed for the square 
forcing case impossible. In the current analysis we will 
therefore content ourselves with scanning the parame- 
ter space numerically to show that, for a wide range of 
systems with triangular forcing, the envelope of the res- 
onance tongues has only one minimum. In order to do 
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this, we first note that the problem now depends on only 
two nondimensional parameters: Ki and K%. We also note 
that we can write the analytic solution of 153|) and use 
(|48|l and l|50[) to obtain a (very long) explicit algebraic ex- 
pression for the envelope of the resonance tongues, which 
we label r' rl (fc/i) but do not reproduce here because of its 
length. By evaluating c^r' rl (fc/i) at 10 3 points between 
k = and k- values that reach an asymptotic regime, us- 
ing approximately 10 4 different (logarithmically spaced) 
combinations of the parameters n\ £ [10 _6 ,10 1 ] and 
k 2 £ [1CT 5 , 10 3 ], we find that r* rl (fc/i) is always a smooth 
function with positive concavity. This strongly suggests 
that Tg rl (fc/i) has only one minimum and that the conjec- 
ture also holds for triangular forcings. 

Finally, for a sinusoidal forcing f(u)t) oc cos(wi) one 
can only express Tk(x) m terms of an integral equation 
which cannot be explicitly solved. The work in |28l |29| . 
however, shows that T s e m (kh) again appears to have only 
one minimum for any combination of parameters. 

The results presented above relate the shape of the 
forcing function to that of the envelope of the resonance 
tongues. In particular, they support the conjecture that 
only a forcing with more than two extrema per cycle can 
generate a tongue envelope that has more than one mini- 
mum. A full proof of this conjecture would be of interest 
not only as a mathematical result, but also as a guide for 
engineering surface patterns. It would imply, for exam- 
ple, that only forcing functions that have this characteris- 
tic can display bicritical points involving non-contiguous 
resonance tongues. 

VI. DISCUSSION AND CONCLUSIONS 

We have presented a new approach for studying the 
effect of the shape of the forcing function on the Faraday 
linear surface wave instabilities. Through a numerical, 
experimental and analytic investigation, we have estab- 
lished a relation between the number of extrema in the 
forcing function and the number of minima that can ap- 
pear in the envelope of the resonance tongues. This ap- 
proach does not rely on a multi-frequency description of 
the forcing function. It therefore allows us to consider 
forcings that cannot be defined by the superposition of 
a few sinusoidal terms, but that can excite surface wave 
instabilities in new ways that could lead to a greater con- 
trol of the surface patterns. 



The analysis that we have carried out provides new in- 
sights for understanding the effects of the energy feeding 
mechanism in pattern forming systems. Indeed, we use 
the lubrication approximation to reduce the system to 
one degree of freedom and then apply the WKB method, 
which neglects the fast oscillations by integrating their 
net effect over the different forcing segments. By doing 
this, we achieve a description that is somehow similar to 
the simple mechanical analogies (with balls, springs and 
pcndula) that are used in reduced dimensionality mod- 
els of parametric resonance. In this context, it would be 
interesting to try to relate the simplified dynamics that 
the WKB calculations furnish for each wavenumber to 
the forcing strength required to reach its corresponding 
instability threshold. Furthermore, it may be possible 
to follow a similar approach to study the effects of the 
forcing mechanism in other fluid regimes or even in a dif- 
ferent system, such as the granular Faraday experiments 
where strongly non-sinusoidal forcings is the norm 

From an analytical perspective, various additional con- 
nections between the forcing shape and the resonance 
tongues could be obtained by developing the implicit re- 
lations established here. We expect to be able to achieve 
this by adequately choosing a reduced set of forcing func- 
tions and using the right approximations. Obtaining 
these additional connections could lead to a better un- 
derstanding of the inverse problem, in which the forcing 
function would be tailored to achieve a given instability. 

From an experimental perspective, the lubrication 
regime in which our analytic results are obtained has not 
yet been widely explored. This is not due to any fun- 
damental limitation but rather to technical difficulties, 
mainly in achieving high enough accelerations at low fre- 
quencies and having a large enough container for the sur- 
face patterns to develop. However, given that we obtain 
good analytical approximations in this regime, we hope 
that new experiments will explore this regime. This, to- 
gether with an extension of our analysis to consider non- 
linear effects, would allow an exploration of the patterns 
that can be formed by the linear instabilities achieved 
through the forcing function control. 
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